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DEGENERATE DIFFUSIONS ARISING FROM GENE 
DUPLICATION MODELS 

By Rick Durrett^'^ and Lea Popovic^'^ 

Cornell University and Concordia University 

We consider two processes that have been used to study gene 
duphcation, Watterson's [Genetics 105 (1983) 745-766] double reces- 
sive null model and Lynch and Force's [Genetics 154 (2000) 459-473] 
subfunctionalization model. Though the state spaces of these diffu- 
sions are two and six-dimensional, respectively, we show in each case 
that the diffusion stays close to a curve. Using ideas of Katzenberger 
[Ann. Probab. 19 (1991) 1587-1628] we show that one-dimensional 
projections converge to diffusion processes, and we obtain asymp- 
totics for the time to loss of one gene copy. As a corollary we find 
that the probability of subfunctionalization decreases exponentially 
fast as the population size increases. This rigorously confirms a result 
Ward and Durrett [Theor. Pop. Biol. 66 (2004) 93-100] found by sim- 
ulation that the likelihood of subfunctionalization for gene duplicates 
decays exponentially fast as the population size increases. 

1. Introduction. Studies have shown that a surprisingly large number of 
duplicated genes are present in all sequenced genomes, revealing that there 
is frequent evolutionary conservation of genes that arise through local events 
that generate tandem duplications, larger-scale events that duplicate chro- 
mosomal regions or entire chromosomes, or genome-wide events that result 
in complete genome duplication (polyploidization). Analyses of the human 
genome by Li (1980) have revealed that at least 15% of human genes are in- 
deed duplicates, with segmental duplications covering 5.2% of the genome; 
see Bailey (2002). For more see the survey articles by Prince and Pickett 
(2002) and Taylor and Raes (2004). 
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Gene duplications are traditionally considered to be a major evolutionary 
source of new protein functions. The conventional view, pioneered by Ohno 
(1970) holds that gene duplication produces two functionally redundant, 
paralogous genes and thereby frees one of them from selective constraints. 
This unconstrained paralogue is then free to accumulate neutral mutations 
that would have been deleterious to a unique gene. The most likely out- 
come of such neutral evolution is for one of the paralogues to fix a null 
mutation and become a pseudogene, or to be deleted from the genome. A 
less frequently expected outcome is the fixation of mutations that lead to 
a new function. However, mutations that lead to such novel gene functions 
are extremely rare; see Walsh (1995), so the classical model predicts that 
few duplicates should be retained in genomes over the long term. 

In Force (1999), introduced a new explanation of the preservation of dupli- 
cates called the duplication-degeneration-complementation (DDC) model. 
To explain, suppose that a gene performs two functions under the control of 
two different regulatory proteins that bind to the DNA upstream from the 
gene. 

fir fir flc 



In the drawing the large rectangle is the gene, typically several hundred 
nucleotides, while the two small rectangles are the transcription factor bind- 
ing sites, typically about 10 nucleotides. It is supposed that mutations which 
cause loss of a regulatory site happen at rate fir while those which cause 
the gene to completely lose function happen at rate fic- In order to have the 
outcome called subfunctionalization in which the two genes specialize to do 
different functions, the first event must be a loss of a regulatory unit. 



S L 



After this occurs, mutations in the indicated regions lead to inactivation 
of one gene copy, /, subfunctionalization, S, or are lethal L since one of the 
functions is missing. It follows that the probability of subfunctionalization 
in a one lineage is 
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If we make the simple assumption that fir = this probabihty is 2/9, but 
if we observe that the gene region may easily be 30 times as large as the 
regulatory elements and set /ic = SOfir then the probability is 1/512. 

Lynch and Force (2000) did simulations to assess the probability of sub- 
functionalization in a population with N diploids. Taking /ic = 10~^, which 
represents a gene with 1000 nucleotides and a per nucleotide per generation 
mutation rate of 10~^ they found that for fir/fJ-c = 3, 1,0.3,0.1 the probabil- 
ity of subfunctionalization, P{S) stayed constant at the value predicted by 
(1) until a population size of about 10,000 and then decreased rapidly to 0. 

In a Research Experiences for Undergraduate project at Cornell, 
Ward and Durrett (2004) investigated the probability of subfunctionaliza- 
tion further using simulation and simple analytical computations. They 
found that if one plots the logarithm of P{S) versus population size N, 
instead of Lynch and Force's plot of P{S) versus log A^, one finds a straight 
line indicating exponential decay of P{S). This conclusion holds for all four 
combinations of haploid and diploid organisms, and linked and unlinked 
duplicates. 

Most multicellular organisms are diploid, but since simulations of 
Ward and Durrett (2004) show that the haploid and diploid situations are 
similar, we will stick with the simpler haploid case here. The linked case, 
where there is no recombination between the two gene copies, is relevant 
to tandem gene duplication in which unequal crossing over between two ho- 
mologous chromosomes leads to a small region of the chromosome present in 
two copies. Ward and Durrett (2004) show that in the haploid linked case, 
the stochastic system converges to a system of differential equations, which 
always ends in a state in which one gene copy is lost. Since, for large N, 
the stochastic system is a small perturbation of this deterministic system, 
it follows from results of Friedlin and Wentzell (1988) that the probability 
of subfunctionalization decreases exponentially in N. 

The unlinked case, in which there is free recombination between the two 
loci, occurs when copies reside on different chromosomes. This case is rel- 
evant to organisms such as the frog Xenopus laevis [Hughes and Hughes 
(1993)], teleost fishes [Bailey, Poulter and Stockwell (1978), Takahata and 
Maruyama (1979), Kimura (1980), Li (1980)] and yeast [Wolfe and Shields 
(1997), Seoighe and Wolfe (1998, 1999)], which have undergone whole genome 
duplication or tetraploidization, followed by a return to diploid inheritance. 
In these species, more than 25% of the duplicated genes have been preserved. 
Using simulation Ward and Durrett (2004) also made the startling discovery 
that in the haploid unlinked model the population frequencies of the various 
genetic states, which should be a six-dimensional diffusion, stayed near a 
one-dimensional curve of equilibria for the deterministic model. 

The main point of this paper is to prove that observation. Our proof will 
show that in directions perpendicular to the curve of equilibria, the system 



4 



R. DURRETT AND L. POPOVIC 



is a small perturbation of a dynamical system with a strong push toward the 
curve. Since the point corresponding to subfunctionalization is not near this 
curve, an appeal to the results of Friedlin and Wentzell (1988) shows that the 
probability of subfunctionalization decays exponentially rapidly as the pop- 
ulation size grows. This rigorously confirms the results Ward and Durrett 
(2004) found by simulation. Our proof gives no insight into the size of the 
constant in the exponential decay. Simulations suggest a value of 0.0005 
for the diploid unlinked case, so for a population size in the millions, the 
probability of subfunctionalization will be extremely small. 

A similar conclusion can be made about models in which the two gene 
copies control more than 2 functions which are regulated by different tran- 
scription factors. In the case of 3 functions the diffusion model becomes 
14-dimensional. The equations for frequencies of all functional factors will 
again produce a one-dimensional curve of equilibria. However, since it be- 
comes unwieldy to compute this curve or show that the 13-dimensional ma- 
trices associated with the linearization have all negative eigenvalues, we shall 
not try to prove this case. 

As is often the case in work of this type, one can be certain of the truth 
of the mathematical fact, but cannot conclusively demonstrate its relevance 
to biology, because the result is based on assumptions which may or may 
not hold. In our case, some suspicious assumptions are: 

1. constant population size, 

2. mutations have no selective advantage, 

3. the duplicated gene is present in all individuals. 

Taking these in the order given, two referees have pointed out that the 
duplications that led to many gene families in Drosophila may have occurred 
well before the divergence of insects and chordates, perhaps in conditions 
of smaller population size. Of course, such gene families must be shared by 
all of these species, so this explanation applies only to the small number 
of families shared by all of these species. A second objection relevant to 
tetraploidization events is that these events initially involved only a single 
individual, so subfunctionalization may have occurred when the population 
is small. A third objection relevant to Drosophila is that they undergo large 
seasonal fluctuations and have spatial structure, so mutations may arise 
when the population size is large and become fixed in a subpopulation when 
it becomes small. Thus subfunctionalization can occur when the population 
size is small but this complements our result that it is unlikely when the 
population size is large. 

As for the second assumption, we have assumed, following the proposers of 
the subfunctionalization explanation, that the mutations involved are neu- 
tral. In some cases, subfunctionalization means that the two genes are ex- 
pressed in different tissues, and further mutations can improve the perfor- 
mance of the proteins. For a concrete example, consider the genes CYPBl 
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and CYPB3 of the Black Swallowtail studied by Wen et al. (2006). These 
genes, which have adapted to metabolize furanocoumarin found in differ- 
ent host plants are unlikely to be ancient duplicates common to all insects 
and chordates, but show a large number of amino acid altering mutations 
indicative of the action of positive selection (see Figure 2 in the paper cited). 

Having admitted that positive selection is possible, we should also add 
that results of Walsh (1995) mentioned earlier show that positive adaption 
explains very few instances of gene duplication. As one referee pointed out, 
the probability of subfunctionalization may indeed be small, but there are 
many genes and many attempts, so there may be many successes anyway. 
One could object to this explanation by noting that if the subfunctional- 
ization probability was 5% then there would be 19 failed duplicates in the 
genome for each success, and this is not observed. However, this objection 
vanishes if the duplications occur in a single individual and only become 
fixed in the population if it provides some benefit, since in this scenario 
failed attempts do not accumulate in the genome. 

In summary, our result does not sound the death knoll for the popular 
concept of subfunctionalization. However, it does indicate that in order for 
it to be at work in a population which currently has a large population size, 
then it must be aided by positive selection or demographic factors. 

2. Gene duplication models. In this section we give precise statements 
on how degenerate diffusions arise in two distinct gene duplication models. 
Sections 3 and 4 contain proofs of Theorems 1 and 2 for Watterson's model 
in 2 dimensions. In Section 5 we give the proof of Theorem 3 for Lynch and 
Force's six-dimensional subfunctionalization model. 

2.1. Watterson's model. To prepare for the complexities of the six-dimen- 
sional model, it is useful to start with Watterson's (1983) double recessive 
null model, which can be thought of as a special case of the subfunctional- 
ization model in which the gene has only one function. Watterson considers 
a Wright -Fisher model with unlinked loci and diploid individuals. Unlinked 
loci means that the two gene copies are free to mutate independently and 
dependence only comes through viability of the gametes. In generation n we 
have 2N letters that are either A (working copy of gene 1) or a (nonfunc- 
tional copy) and 2N letters that are either B (working gene copy of gene 
2) or b (nonfunctional copy). To build up generation n + 1 we repeat the 
following procedure until we have successes: 

• Pick with replacement two copies of gene 1 and two copies of gene 2. 

• An A, that is picked may mutate to a with probability fi. Likewise a B, 
that is picked may mutate to b with probability The reverse mutation, 
which corresponds to undoing a specific deleterious mutation is assumed 
to have probability 0. 
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• If the result of our choices after mutation is not aabb then this is a success, 
and the new individual is added to the collection. 

Note that after forming the new individual we do not keep track of the copies 
that reside in a single individual. This standard practice is referred to as the 
assumption of a "random union of gametes." 

Letting x and y be the frequencies of a and b alleles, Kimura and King 
(1979) derived the diffusion limit of this model to be 

L,f = -x{l -x)g^ + 2iV(l - x)(/. - xV)^ 

o 

To explain the coefficients in (2) we note that if we have a one locus 
Wright-Fisher model in which mutation from A to a occurs at rate ^ and 
allele a has a selective disadvantage s then when time is run at rate 2N the 
diffusion approximation is 

1 df 

In our case an A allele has fitness 1 and a allele has fitness if it is paired 
with another a and two 6's, so the selective disadvantage of an a allele is 
s = xy^. 

This diffusion limit in (2) is unusual since it does not assume that s and 
/i are of order 1/N . Since all individuals have fitness 1 or 0, it is not sensible 
to make this assumption about s, but one can, as Watterson did, assume 
4A^/i 9. By using arguments that were clever but not completely rigorous, 
Watterson (1983) concluded that the mean time until loss of ^ or S had 
mean approximately 

N[\og{2N)-i,[e/2)], 

where "0 is the digamma function. Here we will give a simple proof of a 
version of his result. In our approach, we will assume that // is a constant. 
Apart from simplifying the mathematics, our motivation is that while mul- 
ticellular organisms have different population sizes, most have a mutation 
rate per nucleotide about 10~^, so if we assume genes have 1000 nucleotides, 
/x = 10-5. 

To state our results, we begin by observing that solutions of the ordinary 
differential equation (ODE) 

^ = (l-x)(/i-x^y2). 



(3) * 



^ = (i-j,)(^-iV) 
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Fig. 1 . Picture of the state space for the Watterson model showing the curve of equilibria 
x^y^ = /i and the flow lines for the ordinary differential equation, which are the level curves 
for the projection map. 

have {1 — yt) / {1 — Xt) constant, so the solution moves along the line through 
(1,1) until it hits the equilibrium curve xy = y^. Since the drift in the 
diffusion is 0{N) while its variance is only 0(1), it should not be surprising 
that the diffusion initially remains close to the ODE. 

Theorem 1. Let Zt = (Xt, Yt) be the diffusion in (2) and the solution 
of (3) both starting from {Xq^Yq). Let < e < There is a constant 7, 
which depends on e so that if N is large then for all (Xq, Iq) G [e, 1 — e]^ 



As the reader may have noticed this is not what the frequencies in the 
discrete Wright-Fisher model do on this time scale. They take significant 
jumps on the way to the curve. With more effort we could prove our results 
for the Wright-Fisher model, however, for simplicity, we will, as Watterson 
did, study the diffusion process. 

With Theorem 1 established, it suffices to consider what happens when the 
diffusion starts near the curve. Given {x,y) we define by (1 — y*)/(l — 

X*) = (1 — ?/)/(l — x) and x*y* = yfJL. In words, y) = (x*, y*) is the point 
on the equilibrium curve which we would reach by flowing according to the 
ODE starting from See Figure 1 for a picture. 




sup 

■0<t<7logAf/Af 
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Theorem 2. Consider the diffusion {Xt,Yt) in (2). Let r = inf{t : Xt = 
lorYt = l}he the time to loss of A or B. Let < 6 < 1/2. Suppose \fi — 
^0-^0^1 — ■ Then if N is large, with high probability we have \ix — X^Y^\ < 
2N~^ for all t <t. Also, as N ^ oo the process X^ — Yf* converges in distri- 
bution to a diffusion process on the curve {xy = -^//I}, from which we obtain 
Eqt ~ 2Nc2{fi). 

The coefficients of the hmiting diffusion process can be exphcitly calcu- 
lated by applying Ito's formula to h{^{Xt,Yt)) where h{x,y) = x — y, but 
they have very complicated formulas. Figure 2 shows the drift b{z), variance 
a{z) of the limiting process and the quantity —2b{z)/a{z) which appears in 
derivative of its natural scale 

s'(y)=exp^— y —2b{z)/a{z)dz 

when fi = 10~^. The Green's function is shown in Figure 3. Recalling Eqt = 
2 J^''^ G{0,x)dx and integrating numerically we see that in our concrete 
case C2(/i) = 6.993302, so recalling how time is scaled we have Eqt ~ liN 
generations. 

The result can be explained by noticing that the distance of the diffusion 
from the curve (j){Xt,Yt) = {fj, — XfY^)'^ is a Lyapunov function insuring 
exponential stability for the ODE system (3). Near the curve the stability 
of the drift beats the tendency of the diffusion to deviate from the curve, 
and forces (Xt,Yt) to follow where the flow would take it on the curve. 

Watterson's argument proceeds by making a change of variables p = 
2N^/'^xy (we have s = 1 so his S = N) and 



x-y 
x-y 

\ — X 



li x>y, 
\1 x<y. 



Fig. 2. 



-1.0 -0,5 0.0 0,5 1.0 

'usion coefficients for the limiting diffusion in the Watterson model. 
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Fig. 3. Green's function for the limiting diffusion in the Watterson model. 

Since his /i = 9 /AN, — 9 measure the displacement from the curve x^y^ = 
/i, and — 1 < ?7 < 1 gives the position along the curve. The change of variables 
that defines rj is not twice differentiable when x = y but by smoothing and 
passing to the limit this is not a problem. 

The major unresolved difficulty with Watterson 's argument is his assump- 
tion that p, which runs on a faster time scale than rj, stays in equilibrium 
as rj evolves. A significant technical problem in justifying this is that when 
rj = the diffusion for p does not have a stationary distribution. Indeed, 
it is the time spent with rj near that results in the iVlogA^ asymptotics. 
Of course in his situation p = 0{1/N), so the situation is somewhat more 
complicated. Our proof, in contrast, uses the smooth change of variables 
{x,y) X* — y* , and to obtain our limit all we have to show is that the 
displacement from the curve stays close to 0. 

2.2. Subfunctionalization. We turn now to the subfunctionalization model 
for unlinked loci in haploid organisms. In Watterson's model, the discrete- 
time Wright-Fisher model leads to a process that makes deterministic jumps 
off of the curve of equilibria. To avoid this, we will consider the Moran model 
version in which individuals are replaced one at a time, and in which mu- 
tation happens independent of reproduction. The biology behind the last 
choice is that we are considering the pool of gametes that exist at time t in 
a population with overlapping generations. 

To introduce our model, consider first an infinitely large population for 
which the allele frequency dynamics will be deterministic. Let 3 = 11, 2 = 10, 
1 = 01, and = 00 denote the four possible states of each gene copy, where 
1 and indicate presence or absence of the two functions, and let Xj and yj 
denote the frequencies of states i and j at the first and second copy with 
xq = 1 — X3 — X2 — xi, and yo = 1 — ys — y2 — Hi- To simplify we will assume 
fJ-r = IJ'c = b. Let 

(4) w = X3 + y3- X32/3 + Xiy2 + X2yi 
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be the mean fitness, that is, the probabihty the new individual chosen to 
replace the old one is viable. To explain the formula for w, we note that 
if either gene is in state 3, an event of probability x-^ + — x-^y-^, then the 
offspring is always viable, whereas if neither gene is in state 3, the only viable 
combinations are (1,2) and (2, 1). We are assuming the copies are unlinked 
so the events are independent. 

The diffusion limit of this model is 

1 9^/1 9^/ 



(5) 



2 .,,5:2,3 ^^^^^^ 2^^2,3 ^y^^y^ 

-2NF{xs,X2,xi,y3,y2,yi) • V/, 



where F : 1— > is the vector field determining the evolution of the ODE 
system 

dx^/dt = —X3W + X3 — 36x3, 

dx2/dt = -X2W + X2{y3 + yi) + bx3 - 26X2, 

dxi/dt = -xiw + xi{y3 + 1/2) + bx^ - 2bxi, 

(6) 

dys/dt = -ysw + 2/3 - 361/3, 
dy2/dt = -y2W + y2{x3 + xi) + by-s - 2by2, 
dyi/dt = -yiw + yi{x3 + X2) + by-s - 2byi. 
If we let a = 1 — 36, then the equations for X3 and 2/3 become 



(7) 



dX3 

— = X3{a-w), 
dy3 , ^ 



so the first and fourth equations for an equilibrium reduce to the single 
equation w = a. Thus, if things are nice we will have a one-dimensional 
curve of equilibria. 

To find one set of solutions we can begin by investigating the case in which 
X2 = x\ = X and y2 = Vi = V which gives us four equations 

dx3/dt = —X3W + X3 — 36x3, 
dx/dt = —xw + x(y3 + y) + 6x3 — 26x, 

dy3/dt = -y3W + 2/3 - 36^3, 

dy/dt = -yw + y{x3 + x) + 6^3 - 2by, 

which after some algebra can be explicitly solved. See Figure 4 for a graph 
of the solutions for X3,y3,x,y in the special case 6 = 0.001. 
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It turns out that, in addition to the curve of equihbria, there are no other 
equihbria except for the fixed points of the dynamical system near X2 = 1, 
yi = 1 and xi = 1, y2 = 1, which correspond to subfunctionahzation. 

Equation (7) imphes that if 1/3(0) /x3(0) =r then we will have 1/3 (t) (t) = 
r for all t. Given z = {x3,X2,xi,y3,y2,yi), let z* = be the point on the 
curve of equilibria with y^ /x%^ = 2/3 /x^ . Numerical results show that starting 
from any z the ODE will converge to ^(z), but we do not know how to prove 
this for the dynamical system, so we will only consider the situation when 
the process starts close to the curve. 

Theorem 3. Suppose that the mutation rate h < 0.01. Let r = inf{t : Xo(t) = 
1 or YQ{t) = 1} be the time to loss of gene 1 or gene 2. Suppose \Zq — Zq \ < 
. Then if N is large, with high probability we have \Zt — Zl\< IjN^I^ 
for all t <T. Also, when run at rate 2N the process X^{t) — Y^{t) converges 
in distribution to a diffusion process on the curve of stable equilibria of (6), 
from which we get that Eqt ~ 2Nc2,{b). 

The inspiration for our proof comes from Katzenberger (1991), who found 
conditions for a sequence of processes to converge to a diffusion on a sub- 
manifold of the original state space. While his result gives us the plan, it 
does not help much with the details. The hard work is to obtain enough 
information about the curve of equilibria to be able to check that the lin- 
earization of the drift in the five-dimensional space perpendicular to the 
curve has all eigenvalues with negative real part. (It is at this point we need 
the assumption that b < 0.01.) Linear ODE's of this type have Lyapunov 



b=0,001 



o 




C3 
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0-0 



0,2 



0-4 



06 



0,6 



1,0 



Fig. 4. Picture of the solution of curve of equilibria in the subfunctionalization model 
where X2 = xi = x,y2 = yi = y and the mutation rate has value b = 0.001 — the values of 
{yz,x,y) are shown against xz which is on the horizontal axis. 
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functions, so patching these together and using computations in Section 3 
of Katzenberger (1991) gives a Lyapunov function near the curve. 

Once we know the diffusion stays near the curve, the coefficients of the 
hmiting diffusion can again be computed in terms of the derivatives of 
^{X-s, X2, Xi,Y3,Y2,Yi) and h = — by using Ito's formula. Figure 5 
shows the drift b{z), variance a{z) of the hmit of X^ — Y3* and the quantity 
—2b{z)/a{z) which appears in derivative of the natural scale when h = 10~^. 
Even though the model is quite different, the curves are similar to those in 
Figure 2. The Green's function is shown in Figure 6. Integrating we see that 
in our concrete example 03(6) = 3.284906, so for the process run at rate 1 we 
have Eqt ~ 6.5A^. In the case of Drosophila who have an effective population 
size of = 10^ this is 6.5 million generations or 650,000 years. 

3. Proof of Theorem 1. We will prove this result by using a standard 
estimate from the Picard iteration proof of the existence of solutions of 
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-1.0 -0.5 0.0 0.5 1.0 

Fig. 5. Diffusion coefficients for the limiting diffusion in the suhfunctionalization model. 



0-0 0.2 0.4 06 0.8 1.0 
Fig. 6. Green's function for the limiting diffusion in the suhfunctionalization model. 
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stochastic differential equations (SDE). We first consider a time rescaling 
Zt = {Xt/2N ■,'^t/2N)i so that the coefficients of the SDE for Zt become 

h{z) = ((1 - - xV), (1 - - a^V)), 



yfx{l-x)/2N 



^y{l-y)/2N)- 

Let he{z) and (7e{z) be the formulas that results when x and y are replaced 
by (e/2) Va; A (1 -e/2) and (e/2) V y A (1 - e/2). We do this so that be and 
\/2Nae are Lipschitz continuous with constant K^. 

Let Z° be the solution of the ODE starting from {Xo,Yq) £ [e,l — e]^ 
with e < fi. Since along the boundary of the square, the drift points into the 
square except near the upper left and lower right corners, it is easy to see 
that 

ZO-ZO= fbe{Z's)ds. 

Jo 

To begin the iteration to produce a solution of the modified SDE let 

Z^ = Z^+ fa,iZ^)dB,+ fbe{Z^)ds, 
Jo Jo 

where Bg is a standard two-dimensional Brownian motion. Since u{l — u)< 
1/4 for all G [0, 1], the maximal inequality for submartingales implies 

e( sup \Zl - Z^rt < AE\Z} - Zfp < t/N. 

To continue we will use (2.3) on page 185 of Durrett (1996). 

Lemma 1. If we let 

Ut = Uo+ f a,{Us)dBs+ f\s{Us)ds, 
Jo Jo 

Vt = Vo+ t ae{Vs)dBs+ f\e{Vs)ds, 
Jo Jo 

then for B = {8T + 64)^^ 

(8) sup \Us-VsA <2E\Uo- Vol'' + BE r \Us-Vs\''. 

\0<s<T ) Jo 

Letting Z^ = Z^~^ ^ it follows by induction that for > 1 
e{ sup \Z^,-Zt'?\<^'^~'^' 



0<s<T 



N k\ 
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Using (2.6) on page 188 of Durrett (1996), we see that the solution Zt has 



^fc-i j^fc \ 2 1 ((8T + 64)i^£)* 



The form of the right-hand side derives from the fact that it is the 1? norm 
11^ II 2 = (-EX^)^/^ that satisfies the triangle inequality. The right-hand side 
increases rapidly with T, so we take T to be a constant to end up with a 
bound of the form Ct/N . 

To extend the comparison to longer intervals of time, let Z^ ' and Z^ ' 
be solutions of the ODE and SDE starting from Zj,, and let Z^' be the 
solution of the SDE starting from Zt- It follows from (9) 

sup \Zl'^ - Zl^^\A <Ct/N. 

\0<s<T J 

Using (8) with Gronwall's inequality [see (2.7) on page 188 of Durrett (1996)] 



e( sup \Zl''-Zl'r)< 

\0<s<T / 



2Ct 



N 

Combining the two estimates using the triangle inequality gives 
e( sup \zl'^-Zl^°\A<^{l + V2^f. 



0<s<T 



N 



To continue for A; > 2 let Z^'^ and be solutions of the ODE and SDE 
starting from Z^j, = Zj,~^'^, and let Z^''^ be the solution of the SDE starting 
from Z^^'"^ . We will show by induction that 

(10) e( sup |ZP-Z,^'T) <^fE(2e^^H'. 



0<s<T 



N \^ 



Our previous result shows that this hold for k = \. It follows from (9) 



E[ sup |Z^i-Z^0'2 

vO<s<T 



^ <Ct/N. 



Using (8) with Gronwall's inequality, and our induction assumption we get 



E ( sup I Z 

^0<s<T 



Combining the two estimates gives (10). If we choose T so that 2e^^ 
then 

[fy; ' J - [1 - (2e^^)-V2]2 
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and summing over k= 1, . . . ,m gives 

e( sup \Z.-Z^A<^-^='-^CrN-^/^<N-y^ 

\0<s<mT J iV 1 - 1/4 6 

when we choose m = (log A^)/ (3 log 4). This proves Theorem 1 with 7 = 
T I (6 log 4) and time rescaled back to its original rate. 

4. Proof of Theorem 2. The proof of Theorem 2 is given in four steps. 
We first find the map <I>(rE,y) that takes a point [x^y] to where the ODE 
trajectory would take it on the curve; and apply Ito's formula to find the 
SDE of ^(Xt^Yt). We then use the exponential Lyapunov function for the 
ODE to show that the diffusion stays close to the process and 
finally use its SDE to obtain the limiting process. 

4.1. Projection map of the trajectory. Given x,y in the unit square the 
ODE will take it to the point x* ,y* such that x*y* = while (1 — — 
x) = {1 — y*) / {1 — X*) . Solving for x* gives a quadratic equation 

(1 - y){x*f - (x - y)x* - 07(1 -x) = 

and the root we want is 




Ifx — y Ux — yY 1 



{l-yf ^"^l-y 

To get ready to differentiate this, we note {x — y)/{l — y) = l — {l — x)/{l — y) 
so we can write 



fl-x\ , , (l-M) + w'(l-u)2 + ^6n 
X =9\-^-—^y where = ^ . 

To start to understand the properties of g we note that 

/ (x,y) M = (l-x)/(l-y) 9{u)\ 
(1,^) 1 

1 ^1/4 

V (^,1) 00 ^) 

where for the last evaluation we note that g(u) ~ ^[(1 — u) + ti — (1 — 2y^)] 
for large u. Differentiating gives 
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d'^x* f,fl — x\ 1 



9V_ „ n-x \ {1-xf , n-x \ 2{i-x) 

where 

g'{u) = + i((l - uf + A^mr^/\-2{l - u) + A^)), 



2\ ' 2 

'1/2 



"(.,\ - _ 1 ((1 _ uf + 4V^n)"3/2(-2(l - n) + 4V^)2 



4.2. /to's formula for the projection map. Now consider the behavior of 
= (^[Xt,Yt). We conclude from Ito's formula that 



1 

+ 2 



xy,(i-y,)(is. 



Using (2) to get expressions for dXg^dYs, the drift terms in the first two 
integrals which contain factors of 2N cancel leaving us with 



1 - yj (1 - 

from which we see that the variance process is 



{X*)t-{X*)o= l\'(L^y—±—Xs{l-Xs)ds 







For computing Er, it is convenient to have a symmetric diffusion so we 
will consider X^ — Y^ . Analogously to x* we have 



{y-x) + ^{y-xf + 4^{l-x){l-y)^ n_y 
2(1 -x) ~^\l-x 
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SO 



dy \1 — X 

^ = n' (— y 

dx 

_ 

dy^ 

dh* _ 
dx"^ 



-1 



1 — X 

1-y 



l-x J {l-xY' 



1 

1 



1 



X/ (1 — xY 



l-xj {l-xf 



+ 9' 



l-y\2{l-y) 



1 — X ) {\ — xY 



and we conclude that 
cient 



is a diffusion with two times the drift coeffi- 



2b{Xt,Yt) 



(11) 



l-Xt\Xt{l-Xt 



+ 9 



l-Yt 

l-Xt 
l-Y 



(1 



„/ i-yAy,(i-y,) 



XtfY, 



+ 9 



1 - FA (1 - YtY^t 



1-XtJ (l-XtY 



1 - Xtyji - Xt)Yt 

l-Yt) {l-YtY 
l-Yt\2{l- Yt)Xt 



i-Xt) {i-XtY 



and diffusion coefficient 



a{Xt,Yt) = g' 



(12) 



+ 9' 



l-Xt 
l-Yt 

1 - 



Xt{l-Xt) , {l-XtYYt 



+ 



Yt 



l-Xt 



{l-YtY {i-YtY 

Yt{l-Yt) {l-YtfXt 

{i-XtY ^ {i-XtY 



4.3. Staying near the curve. To study the distance from the curve we wih 
consider (j){Xt,Yt) = {n — XfY^'Y. It turns out that is a Lyapunov function 
for the system (3). In other words, (t){x,y) = {fi — x'^y'^Y ^ Oi 4>{^^y) = iff 
(x, y) is a fixed point of the deterministic system (3), and for all (x, y) in the 
neighborhood of fixed points |/z — x'^y'^\ < the change in the direction 
of the strong drift is 

\/(j).F = -4:xy{y{l - x) + x{l - y))0 < -(3(p, 

where /? = inf^-^j^. |^_3,2j^2|<;v-*{42;y(2/(l - x) + x{l - y))} > 0. 
Using Ito's formula we get 



,2Nl3t 



HXt,Yt) = e'""^' (^2NP^{Xs,Ys) + |^2iV(l - X,)(^ - X^' 
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+ |^2iv(i-y.)(/.-xfy,2) 

and collecting terms due to drift 



For {x,y) within a compact set it is easy to see that for any sequence of 
random times and T > the second integral is less than 

C sup 



0<i<TAr 



N 



Jo 



which converges weakly to 0, likewise the third integral converges weakly to 
as well. 

Let = inf{t > : |^ - XfY^^\ > N-^}, then for all < t < the first 
integral is negative, and since (/> > this implies (p{X^/^^N ,Y^/^^n) ^ 0. In 
other words 

i^tArN ,Yt^^N) - {X*^^N ,Y*^^n) =^ 0. 

4.4. Limit theorem. Our final step is to argue that as N ^ oo, X^ — Y^* 
converges to the diffusion with coefficients b{Xf,Y^) and a{X*,Y^*). We 
have shown that for each the drift and diffusion coefficients are given 
by b{Xt,Yt) and a{Xt,Yt) as defined in (11) and (12). Since the coefficients 
b and a are bounded in a neighborhood of the curve of equilibria the se- 
quence of processes is tight, and the weak convergence {X^/^^n jY^/^^n) — 
{X*^^j^ ,Y*^^tf) =^ that we have just shown together with Theorem 5.4 in 
Kurtz and Protter (1991) implies for all < t < we have convergence of 
X^ — Y^ to the desired diffusion process. 



5. Proof of Theorem 3. We follow the same general outline as the proof 
of Theorem 2, but there are some new steps which were trivial in the pre- 
vious proof. First, finding a usable expression for the curve of equilibria 
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requires more work, and to calculate the flow map requires the right change 
of variables to describe the surface along which it is constant. Next, we were 
not able to find a global Lyapunov function in this case, hence we linearize 
the system around the points on the curve, and investigate the behavior of 
the eigenvalues of the linearized system. As the last step we use the con- 
struction of Katzenberger to glue the local Lyapunov functions into a global 
one for the whole neighborhood of the curve. 

5.1. The curve of equilibria. Consider the equations in (6) restricted to 
Xi=X2= X, yi = y2 = y. 

(13) ^^ = -X3W + X3-3bx3 = X3{a-w), 

dx 

(14) — = -xw + x{y3 + y) + 6x3 - 2bx = x{y3 + y -2b-w) + 6x3, 



(16) — = -yw + y(x3 + x) + bys - 2by = y{x3 + x - 2b - w) + by^, 

where w = X3 + ys — x^y^ + 2xy and a = 1 — 3&. Letting [) = a + 2b = l — 
fixed points satisfy three equations 



(15) 




dy 




2^3 + 2/3 - 2:32/3 + 2xy = a, 
x{y3- (3) +xy + 6x3 = 0, 
y{x3 - f3) + xy + by3 = 0. 



(18) 
(19) 



To solve it is convenient to let 



(20) 



r = X3 + 2/3 - X3y3 - a 



2xy. 



We can solve (17), (18), (19) for x and y 




y = 



r - 257/3 
2(x3-/3)' 



Using these in (20) we have 



(22) 



2r(x3 - /3)(2/3 - /3) + (r - 26x3) (F - 262/3) = 0. 



Filling in the definition of F 

2(-X3y3 + X3 + ys - a)(x3?/3 - [3x3 - Pys + /3^) 

+ {-X3y3 + X3(l - 26) + 2/3 - a)(-X3y3 + X3 + 2/3(1 - 26) - a) = 
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and collecting terms gives 

xlvl [-2 + 1] + (x^yg + x^yl) [2/3 + 26] 

+ X3y3[-2/32 - 4/3 + (1 - 2bf + 1] + {xl + yi)[-2/3 + (1 - 26)] 

+ (xs + yg) [2/3^ + 2a/3 - 2a(l - 6)] + [a^ - 2a(3^] = 0, 

which has the form j Cijx\y-'^, where Cjj- = Cj^i. Using /3 = 1 — 6 and a = 
1-36 

C2,2 = -1, 

ci,2 = C2,i = 2/3 + 26 = 2, 

co,2 = C2,o = -2/3 + (l-26) = -l, 

ci 1 = -2(1 - 26 + 6^) - 4(1 - 6) + (1 - 46 + 46^) + 1 = -4 + 46 + 26^, 
co,i = ci,o = 2/32 = 2 -46 + 262, 

Co = - 2a/32 = a(a - 2/3^) = (1 - 36)(-l + 6 - 26^) 

= -1 + 46- 56^ + 66^ 

To solve for 2/3 note our equation has the form do + diys + ^2^/3 = where 
dj = CijX^^, so 



-di ± i/df - 44^2 

<^'' ^5 ■ 

To see which root we want, note that if 3:3 = and x = then w = = a, 
so we want 



-co,i + Jcg 1 - 4co,oCo,2 

ys = 

C0,2 

-2/32 + V4/34 -4(a^ -2a/3^)(-l) 



a = 1 — 36. 



5.2. Projection map. It follows from the equations for the ODE that 

d_y3 ^ ^^dX3 _^ 

2:3 X3 2/3 dt 

so ys/x^ is constant along solutions. This means that if 2/3/2:3 = r then the 
trajectory will flow the point on the curve of equilibria where n = X3 has 

{ru)'^d2{u) + {ru)di{u) + do{u) = 0. 

Differentiating with respect to r, and letting = (u^di)' , where ' indicates 
the derivative with respect to u we have 

2 r 2 du , du du 

2ru d2+r e2— + udi + rei— + cq— = {). 

dr dr dr 
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Rearranging gives 

du 2rv?d2 + udi 

(24) 



dr r'^e2 + rei + cq 
From this it follows that 

d du d I du 

-i—u{y3 X3) = 2 3-3— ^(ys/a^s) = —-f 

dX3 X3 dr ay^ X3 dr 



d^ ^ I ^ _ yl d^u ^2y2,du d^ ^ ^ ^ _ 1 d^u 
(ixg X3 dr"^ X3 dr dy^ Xg dr^ 

To now compute d^u/dr^, we note that du/dr = — f (r) / g{r) with f'{r) = 
2u^d2 + (2re2 + ei)du/dr and g'{r) = 2re2 + ei + (r^e'a + re'i + e'Q){du/dr) 
where is the derivative of with respect to r. 

Let V = y^ and g = X3/7/3 = 1/r. To compute the derivatives of v we note 
that if dj{z) = ^iCijZ^^ then repeating the previous steps leads to 

dv _ 2qv'^d2{v) +vdi{v) 

dq q'^e2iv) + qei{v) + eo{v) 

and it follows that 

d 1 dv d x^dv 

-l—v[x-iyi) = — — , —v{X3y3) = ^ — , 

dx3 ys dq dy^ y| dq 

d^ , . . 1 d^v d"^ . , . Xn d^v 2x3 dv 

-T^v[x3/y3) = — — , -riv{x3/y3) = —— + — — . 
dxi yi dq^ dy^ y| dq^ dq 

To compute d'^v/dq'^, we note that dv/dq = — f2{<l) / 92{<l) with f2{<l) = 2v'^d2 + 
{2qe2 + ei) dv/dq and g2{q) = 2qe2 + ei + {q'^e'2 + qe[ + e'Q){dv/dq). 

5.3. Ito's formula for the projection map. Writing Rg = ¥3(3) / X3{s) we 
conclude from Ito's formula that 

xut)-xm 



u"iRs)-^^ + u'{Rs) 



X3is){l-Xs{s))ds 



Xlis) ''-'^ Xlis) 

The drift terms in the first two integrals cancel leaving us with 
j\'{Rs) • :^^^Xs{s){l-Xs{s))dBl 
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t 



u'iRs) ■ ■j^^^jYsis)il-Y3is))dBl 



so the variance process 

(^3)t-(^3)0 



X3is)il-X3{s))ds 



+ 



u'iRs 



1 



;Y3{s){l-Y3{s))ds. 



Let Qs = X3{s)/Y^{s). Combining these with the corresponding formulas for 
1^3* (t) we see that (t) — Y^ (t) is a diffusion with two times the drift equal 
to 



2b{Xs{t),Ys{t)) 



xm 



XI (t) 



u"{Rt 



xm 



Y3{t){l-Ysit))-v"{Qt 



X3(t)(l-X3(t)) 

1 



Yiit) 
Ysim-Ysit)) 



x,m-x,it)) 



and variance 

a{X3{t),Ys{t)) 

= u{Rtf- 
+ v{Qtf 



^^X,{t){l-X,{t)) + 



1 



,X3(t)(l-X3(t)) + 



x,{t)' 

xm 



Y^{t){l-Y:,{t)) 



Yiit) 



Y^{t){l-Y^{t)) 



5.4. Linearization. In order to show that the diffusion stays close to the 
curve of equilibria, we will investigate the linearization of the ODE around 
these fixed points. Recall from (6) the equations are 

dxs/dt = —X3W + X3 — 36x3, 

dx2/dt = -X2W + X2{y3 + yi) + bx3 - 26X2, 

dxi/dt = -xiw + xi(y3 + 2/2) + 6x3 - 2bxi, 

dys/dt = -ysw + 2/3 - 35y3, 

dy2/dt = -y2W + y2{x-i + xi) + 67/3 - 2hy2, 

dyi/dt = -yiw + yi(x3 + X2) + hy^. - 2byi 
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with w = X3 + y3- x^y-s + X2yi + xiy2. 

We write the coordinates z = {zi,Z2,Z3,Z4,Z5,zq) = {x3,y3,X2,xi,y2,yi). 
If we have a dynamical system of the form dzi/dt = Fi{z) and we let z = 
z* + ep where z* is a fixed point, then comparing terms of order e we see 
that the linearization at the fixed point is 

dpi ^ dFi 

To compute the linearization it is useful to note that Vu; = (1 — y3,l — 
X3,yi,y2,xi,X2) and that at z* we have w = a = l — 3b 

VFi{z*) = -X3VW, 
VF2{z*) = -y3Vw, 

VFs{z*) = -X2VW + (6, X2, ys + yi - 1 + ^ 0, 0, X2), 

VFa{z*) = -xiVw + (6,xi,0,y3 + y2 - 1 + 6,xi,0), 

VF^{z*) = -y2Vw + {y2,b, 0,y2,X3 + xi - 1 + b,0), 

VFq{z*) = -yiVw + {yub,y2,0,0,X3 + X2 -1 + 6). 

On the curve of equilibria X2 = xi = x and y2 = Vi = y, so the four- 
dimensional space p3 =P4, P5 =P6 is invariant for the linearization, as is the 
orthogonal two-dimensional subspace of vectors of the form: (0, 0, u, —u, v, —v) 
Thus our six-dimensional problem decomposes into a two-dimensional one 
and a four-dimensional one. Taking the smaller one first, multiplying the 
linearization on the right the linear map in this subspace is 

(u, v) (ii(y3 + y — 2b — w) — vx, —uy + f (xs + x — 2b — w)). 

To simplify we note that in equilibrium y^+y — 2b — w = -bx^/x and X3 + 
X — 2b — w = —by^/y so the matrix is 

-bx^/x —X \ 

-y -byz/yj' 

For a 2 X 2 matrix M, trace(M) = Ai + A2 and det(M) = AjA2, so the 
real part 3ft(Aj) < if and only if trace(M) < and det(M) > 0. Clearly 
trace(M) = —bx^/x — 6y3/y < 0. The sign of det(A/) = b^x^y^/xy — xy is 
not so clear but as it turns out (for the proof see Appendix). 

Lemma 2. bx^/x > y and by^/y > x so b'^x-^y^/xy > xy. 

Turning to the four-dimensional subspace, we use coordinates (zi, Z2, 23, 24) 
(x3, ?/3, X, y), and the linearization reduces to 

/ -xsil-ys) -X3(l-X3) -2x3y -2x3X ^^ 

jr^ -yall-ys) -y3(l-a:3) -2y3y -2y3X 

— x(l — y3) + 6 — x(l — X3)+x —2xy — bx3/x — 2x^ + x 
V-y(l-y3) + y -y{'^- ^s) +b -2y^ + y -2xy -by^/yj 
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where we have again used y^ + y — 2b — w = — 6x3/2; and X2, + x — 2h — w = 

-bys/y- 

To simphfy, we will change to a new basis vi, V2, v-s, V4. To see how this 
will affect the linearization, note that if V is the matrix with ith column Vi 
then a vector with coordinates q in terms of the new basis, is in the original 
coordinate system 



Pj 



To return to the q system q ■ 



k 

V~^p, so 



Vq. 



dt 



^TVq. 



To cancel the multiples of Vtt; that appear in all of the rows while preserving 
the symmetry in the problem, we choose 

/ l/2;3 -l/?/3 
1/2:3 1/^3 
-x/x^ 

-y/yz 



V- 



\ 







1/ 



V-- 



( X3/2 

-y3/2 

xl2 

\-y/2 



X3/2 
y3/2 

x/2 

y/2 





1 







1/ 



It is easy to check that V 
I 



= 1. The left multiplication simplifies 



-2(1 -y3) 



V 



h 
y 





-2(1 -X3) 

X 

b 




-iy 

-bx-^/x 

y 



\ 





x 

-bys/yJ 



while the right multiplication make things a little messier, 





-2:3(1 -y3)+ 1/3(1 - 
-x{yi + y)/2 
y{x3 + x)/2 



X3 



- 2:3(1 





■ ys) - ysi^ - xs) 

x{y:i+y)/2 
y{x3+x)/2 



- Axy 





-iy 
-bx^/x 

y 





-4a- 

X 

-bys/y, 



This matrix has one zero eigenvector corresponding to the direction along 
the curve. To find conditions which guarantee that all of the eigenvalues of 
the three-dimensional operator on the perpendicular subspace have nega- 
tive real parts, we turn to the Appendix 2 of Murray (1989). Given the 
characteristic polynomial 

P(A) = A" + aiA"-^ + --- + a„ = 

the Routh Hurwitz conditions in terms of the a's are necessary and sufficient. 
For a 3 X 3 matrix M = [niij] these are 

ai > 0, as> 0, 



aia2 - 03 > 0, 
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T - 



tog{-det(M)) 




0-0 0,2 0.4 0,6 0,8 1.0 

Fig. 7. Picture of the of coefficients determining the behavior of the linearization of the 
deterministic drift at the curve of equilibria in the subfunctionalization model. 

which in terms of bi = trace(M), 63 = det(M), and 

62 = ?7l22"^33 + W-limss + miim22 - m23"i32 - "ll2"T'21 " ?«13"l31 

= (miim22 - mi2m2i) + (m22m33 - 771,23^32) + (miim33 - mi3m3i) 

are 

(26) trace(M) < 0, det(M) < 0, det(M) > trace(M)62. 

Note that for this to be possible we must have 62 > 0. 

Before trying to prove these three conditions, we computed them numer- 
ically for the concrete example b = 0.001. In Figure 7 we have plotted the 
values. Since the determinant turns out to be small, we have plotted the 
logarithm of — trace(M), — det(M) and — trace(M)62- 

The eigenvalues of our 4x4 matrix V~^J-V are a and the eigenvalues 
of the 3x3 matrix 




-Ay 
bx-^/x 

y 




It is clear that trace(M) < 0. The determinant is 



det(M) = (-X3(l - y3) - y3(l - x-i) - Axy) 



-xy 



xy 



xjy-i + y) f byz 
2 V y 
xy{y3 + y) , / 




bx3 \ y{xs + x) 
X J 2 



xy{x'i + x) 
2 



so Lemma 2 implies det(M) < 0. 
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For the third condition, since niu < and since Lemma 2 imphes that 

b^xsys/xy > xy, we have 

bx3 fb'^x^y^ \ 

02 = rnn + ^xy{y3 + y)+ I xy 

X \ xy J 

(28) 

mil + '^xy{x3 + x) > 0, 

y 

which is necessary for det(M) > 62trace(M). Once 
Lemma 3. det(M) > 62trace(M). 

Lemma 3 is proved (see the Appendix), we will have verified the third 
condition in (26) and shown that all the eigenvalues have negative real part. 

5.5. Staying near the curve. In order to construct a local Lyapunov func- 
tion in a neighborhood of the whole curve of equilibria with exponential 
convergence properties we follow closely the arguments of Katzenberger. 
We present the outline of the proof and refer the reader to Section 3 of 
Katzenberger (1991) for more details. 

To begin we consider a linear dynamical system 

(29) ^ = ^"(*)- 

If things are nice, for example, if all eigenvalues are distinct, when we let V 
be the matrix with columns equal to the eigenvectors of A, we will have 

A = VAV~\ 

where ^ is a diagonal matrix with An = Aj the eigenvalues of A. Let e*"^ = 
EkLo The solution of (29) can be written 

x{t) = e*^2;(0) = Ve^^V^^x{0). 

Using a superscript T for transpose, a little linear algebra shows 

\\V-'x{t)f = x(0)^(y-i)^e2^V-ix(0) = ^ e'^^\V-'x{0)f, 



11 ) 



so if all eigenvalues have negative real part ||y~^x(t)|| decreases exponen- 
tially. 

To do this in general for a nonlinear dynamical system 

consider the linearized system with A = [ff^(-z*)]ij' along the curve of fixed 
points z* . If all of the eigenvalues of A have negative real part, then there is 
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a > so that they all lie in D{6) = {z G C : |(5z + 1| < 1}. B = 1 + 6A has 
spectral radius p{B) < 1. p{B) = infggQ^^) \Q~^BQ\ so there is an invertible 
matrix V with 



and a little linear algebra shows that decreases exponentially at 

rate a = 1 — p{B). 

We now construct a Lyapunov function in the neighborhood of an arbi- 
trary fixed point x on the curve. Consider first the case that the curve of 
fixed points in the neighborhood W of any fixed point x is of form N r\W 
for some linear subspace A^. Then N = 'Ker^A) and P = Ran(74) decompose 
the whole space, and the spectral radius of B = I + A restricted to P is 
p{B)\p < 1. There is an invertible matrix V on P with ||y~^i?|py|| < 1. 
Extend V to the whole space by setting V = VHp + IIn with projection 
maps Hp, and onto P and A^, respectively. Then ||y~"'^npx(t)|| has an 



Let 9{y) = F{y) - Ay, then De{x) = = 0, and for e > there is 

a neighborhood W CW of x such that ||y"^D0(y)F|| < e for all y £ W, 
and consequently < e||y~^np2;||. Take e < a, and let r = inf{t > 

0:x{t) ^ W'}. Integration by parts on e~^^x{t) gives 



V-'BV\\ < 1. 



Since A = (B — 1)/5 and / and B commute, we have 





hence 



V~'^Upx{t)\\ <e 




1 



npx(s)||+ / e 



('-^)\\V~'^Upe{x{r))\\dr 



J s 




Gronwall's inequality implies 




V-^UNx{t) -V-^Unx{0)\\ 



(31) 





y-inpx(o)||(i-e-("-^)*). 



a — 
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Adding e/(a — e) times (30) to (31) 

\\v-^UNx{t)\\ + ^—\\v-^Upxit)\\ < \\v-'^Unx{o)\\ + ^— ||y-^npx(o)||. 

a — e a — e 

So, if p>0 is such that W" = {x:\\V-^Unx\\ + ^\\V-'^Upx\\ <p}cW\ 
then for any starting point in W" the exit time r = oo. Let 

„, . \\V~^Unx\\ e \\V-'Upx\\ 

f{x) = , g{x) 



p a — e p 

v{x) = g{x) + \J /^(x) + g'^{x), u{x) = (1 — v{x))^ 

and W" = {xe W" : u{x) > 0}. On W" let ^(x) = g^{x)/u{x). Then, 

v{x{t))<v{x{0)), (/)(x(t)) <e-2(°-^)V(x(0)), V<^-F< -2(a-e)(/). 

Note that (j) satisfies > (j){x) = if x is a fixed point, and the reason for 
using the function u{x) in the construction is that 4>{x) = oo for x outside 
W". 

In general let N be the tangent space of the curve of fixed points at x. 
Then = Ker(A) and P = Ran(A) decompose the full space and there is a 
smooth coordinate function r] that will map any x € N to r]{x) G P in such a 
way that x + ?/(x) is a fixed point. Using this coordinate transformation we 
can construct a Lyapunov function in a neighborhood of x in the analogous 
fashion to the one above. 

We can now patch these local Lyapunov functions together over the whole 
curve of equilibria. If (pi is a Lyapunov function constructed for the neigh- 
borhood Wi and (j)2 for W2, define a function onW = Wi U W2 by cp^x) = 
if X is a fixed point and (j){x) = 1/(01 (x) + 4'2{x)) for any other point of W. 
It is easy to verify (p has the same properties as the functions (pi. 

Given a local exponential Lyapunov function (p the proof that the diffusion 
Z{t) = {X3{t),Y3{t),X2it),Xi{t),Y2it),Yi{t)) converges weakly to the curve 
of equilibria and the proof that the limiting process is characterized by an 
SDE derived for ^{Zt) follow by the same arguments as those described in 
Sections 3.3. and 3.4. 

APPENDIX 

Proof of Lemma 2. In equilibrium w = a = 1 — 36, so by (13) 

< bx-s/x = l- b-y3-y, 0< bys/y = l-b- X3- x. 

Thus if can show ya + 2/2 + yi < 1 — 6 (and hence by symmetry X3 + X2 + xi < 
1 — 6) the desired result will follow. To do this we note that w = a = 1 — 3b 



so 



^ = 6y3 + 26y2 + 26yi + x^yo - yo(l - 36). 
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Thus in equilibrium 6(1 — i/q) + b{yi + 1/2) = yo{l — 36 — X3) and it follows 
that yo > V(l - 2^* - a;3) > 6. □ 

Proof of Lemma 3. To check the condition det(M) > 62trace(M) we 
rearrange the formula (28) for the determinant of M = [rriij] (27) 



The trace of M, consists of three negative terms ma, while 62 consists of 
five positive terms which we will refer to as 62^ where j is the order in 
(28). To prove the desired inequality it is enough to find terms smaller (i.e., 
more negative) than the five parts of det(M), di, . . . ,d^, within the fifteen 
products in 62trace(M). The first three are easy: di = 77111623, (^2 = H33622, 
d^ = 77122^25- For the final two we will use the next lemma. 

Lemma A.l. If b < 0.01 then in an equilibrium we have 2/3(1 — X3) > 
x>y whenever X3 > 7/3. 

Using the result of Lemma A.l we see that when X3 > y^, 771-11(622 + 
^25 ) < — ( 1 — 2:3 ) (^22 + ^25 ) < ^4 + c?5 • Symmetry implies that when X3 < 7/3 , 
"111(622 + ^25) < -2:3(1 - y3)(622 + ^25) <di + d5 as weh. □ 

Proof of Lemma A.l. We want to show that if 6 < 0.01 then in equi- 
librium, 



whenever 2:3 < 7/3. Recall that if /3 = 1 — 6 and a = 1 — 36, the fixed points 
satisfy three equations 



For the second inequality we note that taking the difference (34) minus (35): 

x{y3 - P)- y{x3 -P) = Kvs - X3)- 

Since < y^ < f3 = 1 — b (recall that in the proof of Lemma 2 we showed 
7/0 > b) it follows that 



det(M) =mii 




xy \ y X 

y ■ 2xy{x3 + x) - x ■ 2xy{y-i + y). 



(32) 



^3(1 -a;3) >x>y, 



(33) 
(34) 
(35) 



X3 + y3- 2:32/3 + 2xy = a, 
x{y3- P) +xy + 6x3 = 0, 
y{x3 - P) + xy + by3 = 0. 
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To prove the first inequahty in (32) we will show that 2/3(1 — X3) — x is a 
decreasing function h{x3) plus an error of 0(6), and the value of h at the 
point where X3 = 2/3 is positive and 0{Vb). Our equation (23) for 2/3 in terms 
of X3 has the form 



-di + Jdj -4doc?2 

= ' 

where dj = J2i ^"^^ matrix form 

-1 + 46-562 + 63 2 -46 + 26^ 
2 - 46 + 262 _4 + 45 + 262 
-1 2 
with the rows and columns in the order 0, 1, 2. 

d2 = (-l + 2x3-x2) = -(l-X3)2, 

di = {2- 4x3 + 2x1) + + 4x3) + 26^(1 + X3) 

= 2(1 - X3)' - 46(1 - X3) + 26^(1 + X3), 
do = (-1 + 2x3 - xl) + 6(4 - 4x3) - 6^(5 - 2x3) + 6^ 
= -(1 - X3)2 + 46(1 - X3) - 62(5 - 2x3) + 6^. 
If 6 = then dQ + d2 = —di so 

(1 - 2/3) (do - d2y3) =do + diys + ^2^1 
and 2/3 = 1- One can also see this from 

dl - Mod2 = dl + 4(di + d2)a!2 = (di + 2^2)2 = 0. 
Using the definitions of the di, we have 

dl = 4(1 - X3)^ - 166(1 - X3)^ + 1662(1 - X3)^ + 862(1 - X3)^(l + X3) 

- 166^(1 - X3)(l + X3) + 46^(1 + X3)2, 
-4d2do = -4(1 - X3)^ + 166(1 - X3)^ - 462(1 - X3)^(5 - 2x3) + 46^(1 - X3)^ 
Adding the last two equations gives 
dj - Adod2 = 1662(1 - X3)2 + 862(1 - X3)^(l + X3) - 462(1 - X3)^(5 - 2x3) 
- 166^(1 - X3)(l + X3) + 46^(1 - X3)(l - X3) + 46^(1 + X3)2 
= 462(1 - x^f[l + 4x3] - 46^(1 - X3)[3 + 5x3] + 46^(1 + X3)' 
= 462(1 -X3)^[l + 4x3] -M, 
where the mess 

b 3 + 5X3 , 62 (1+X3)2 



(I-X3) (1+4x3) (1-X3)2 (1 + 4x3)' 
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Using 2d2yz = + y(^'i~ ^dod2, we see that 

-2(1 - xsfys = -2(1 - xsf + 46(1 - xg) - 26^(1 + xg) 
+ 26(1 - X3)Vl + 4x3M^/2 

and we have 

ya = 1 - + VI + 4x3) 

1 - X3 



+ n — m ( 1 + 2:3 - (1 - 2:3)71 + 4x3 



(I-X3 
1 - g(x3)6 + r(x3,6)6^. 



Note that g(x3) — > 3 as X3 — > 0, which is consistent with our earher calcula- 
tion that y3 — > 1 — 36 and we have q{x3) > 3 for X3 G [0, 1]. 
To find X using (21), we first compute 

r = y3 + 2:3(1 -2/3) -(1-36) 

= -qb + r6^ + X3(q'6 - r6^) + 36 

= 6(3-g(l-X3)) + r62(l-x3) 

and then 

r - 26x3 6(3 - q{l - X3)) + r62(l - X3) - 26x3 



2(2/3-/3) 2[l-(?6 + r62-(l-6)] 

Multiplying top and bottom by —1/6 

g(l-X3) + (2x3-3)+6r(l-X3) 
^ = 2(g-l)-2r6 ' 

To check the condition 2/3(1 — X3) > x now we note that ignoring terms of 
order 6: 

q{l - X3) + (2X3 - 3) 



y3(i -2:3) - X = 1 - X3 



2(g-l) 



2 + VI + 4x3 + 2x3 - 3 ,^ , 
1 — X3 , • 1 — X3 

2(i+x3 + vmii) ^ 



2(1+X3 + V1 + 4X3) 

It is not immediately obvious that h is decreasing but writing h = (1 
X3)f /{2g), where f,g>0 on [0, 1] we have 

dh f (I-X3) (4/2)(l + 4x3)-^/2 (1 , (4/2)(i I 4a;,)-i/2w 
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The numerator of the final fraction is 

2(1 + 4x3)-^/^ [1 + X3 + VI + 4rc3] - [l + 2{l+4xsy^/\3 + VTT^) 



= (2 + 22:3 - 6)(1 + 4x3)"^/^ + 2 - 2 - 3 - VI + 4x3 < 0, 

so dh/dx3 < 0. To complete the proof now, we will first approximately eval- 
uate h{x3) at the point where X3 = 7/3, and then investigate the errors in our 
approximations . 

Solution when X3 = 1/3. We proceed by guessing and then verifying our 
answer. From the previous calculation 7/3 = 1 — Qb/{1 — x^), where Q ~ 
2 + V5. If X3 = 1 — uVb = ys, then 

2 + V5 



and hence n = v 2 + Vs. 
u 

From (36), since q is large, X3 1, and q = uj^fh 

_ g(l-a;3)-l _ 2 + V5-1 
2g ^ 2-U/V6 ' 

so we have x = y = v\fh with v = -"-^^ . 

2V2+v^ 

To derive this from the equations for the equilibria, we note that when 
3^3 = 2/3 and x = y, (33) implies: 

l-il-xsf + 2x^ = 1 -3b, 

-u^b + 2vH = -3h 

or — 2f ^ = 3. Using (34) and dropping terms of order 6^/^ 

x{x3 - (1 - 6)) + x^ + bx3 = 0, 
-t;n6 + f ^6 + 6 = 0. 
The solution to the quadratic v'^ — uv + 1 = is 

u± \Ju^ - 4 



Taking u = y2 + ^/b and choosing the minus root gives 



2 + V5- V V5 -2 



2 + V5- V(V5-2)(2 + V5) 1 + V5 



2V2 + V5 2V2 + V5 
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For the other equation we note that 

4(2 + ^/5) 2 + 

At X3 = 2/3 we have 

7/3(1 — X3) — a; PS nVo vo = — V > 0. 

2m 2m 

Error analysis. To evahiate the error, we write q = Q/(l — X3) and r = 
i?/(l — 2:3)^ . Taking the difference between the exact expression in (36) and 
the approximation obtained by setting 6 = gives 

^ g(l - X3) + (2x3 - 3) + 6r(l - X3) _ g(l - X3) + (2x3 - 3) 
2(g-l)-2r6 2(g-l) 

_ br[-{l - X3)2{q - 1) + q{l - X3) + (2x3 - 3)] 
4{q - 1 - rb){q - 1) 

bR[-Q-l] 



AiQ-l-Rb/{l-X3)){Q-l) 



We have Q + 1 < 3 + \/5. To bound we begin by noting that if 6 < 0.01 
and 1 — X3 > 6^/^ then 

1 - AfV2 

1 > M> 1 -36/(1 -X3) and hence 0<(1-X3) < 3, 

which imphes < i? < 2 + 3^5 = 8.7082. Since UVb < 1 it follows that when 
1 — X3> 6^/'^ the absolute value of the error is 

^ ,(2 + 3V5)(3 + V5) ^ 19 + 7^5 ^ ^^33^^^ 



To see if this is small enough, we note that m = -y/ 2 + ^5 = 2.05817 so (m — 
1)/2m = 1.27202, and we do have 1.272\/6 > 3.9066 when 6 < 0.01. The last 
detail to check is that when X3 = 6^/^ and 6 < 0.01 

y3 < 1 - 36^/2 + 8.70826 < 1 - 26^/^ 

so the assumption 1 — X3 > 6^/^ we used to bound M is valid. □ 

REFERENCES 

Bailey, G. S., Poulter, R. T. M. and Stockwell, P. A. (1978). Gene duplication in 
tetraploid fish: Model for gene silencing at unlinked duplicated loci. Proc. Natl. Acad. 
Sci. USA 73 5575-5579. 

Bailey, J. A. (2002). Recent segmental duplications in the human genome. Science 297 
1003-1007. 



34 



R. DURRETT AND L. POPOVIC 



DuRRETT, R. (1996). Stochastic Calculus. Probability and Stochastics Series. CRC Press, 
Boca Raton, FL. A practical introduction. MRMR1398879 

Force, A. (1999). Preservation of duplicate genes by complementary degenerative muta- 
tions. Genetics 151 1531-1545. 

Friedlin, M. I. and Wentzell, A. D. (1988). Random Perturbations of Dynamical 
Systems. Springer, New York. 

Hughes, M. K. and Hughes, A. L. (1993). Evolution of duplicate genes in the tetraploid 
animal Xenopus laevis. Mol. Biol. Evol. 10 1360-1369. 

Katzenberger, G. S. (1991). Solutions of a stochastic differential equation forced onto 
a manifold by a large drift. Ann. Probab. 19 1587-1628. 

KiMURA, M. (1980). Average time until fixation of a mutant allele in a finite population un- 
der continued mutation pressure: Studies by analytical, numerical, and pseudo-sampling 
methods. Proc. Natl. Acad. Set. USA 77 522-526. 

KiMURA, M. and King, J. L. (1979). Fixation of a deleterious allele at one of two "dupli- 
cate" loci by mutation pressure and random drift. Proc. Natl. Acad. Sci. 76 2858-2861. 

Kurtz, T. G. and Protter, P. E. (1991). Weak limit theorems for stochastic integrals 
and stochastic diff'erential equations. Ann. Probab. 19 1035-1070. 

Li, W. H. (1980). Rate of gene silencing at duplicate loci: A theoretical study and inter- 
pretation of data from tetraploid fishes. Genetics 95 237-258. 

Lynch, M. and Force, A. (2000). The probability of duplicate gene preservation by 
subfunctionalization. Genetics 154 459-473. 

Murray, J. D. (1989). Mathematical Biology. Springer, New York. 

Ohno, S. (1970). Evolution by Gene Duplication. Springer, Berlin. 

Prince, V. E. and Pickett, F. B. (2002). Splitting pairs: The diverging fates of dupli- 
cated genes. Nat. Rev. Genet. 3 827-837. 

Seoighe, C. and Wolfe, K. H. (1998). Extent of genome rearrangement after genome 
duphcation in yeast. Proc. Natl. Acad. Sci. USA 95 4447-4452. 

Seoighe, C. and Wolfe, K. H. (1999). Updated map of duplicated regions in the yeast 
genome. Genetics 238 253-261. 

Takahata, N. and Maruyama, T. (1979). Polymorphism and loss of duplicate gene 
expression: A theoretical study with tetraploid fish. Proc. Natl. Acad. Sci. USA 76 
4521-4525. 

Taylor, J. S. and Raes, J. (2004). Duplication and divergence: The evolution of new 

genes and old ideas. Ann. Rev. Genet. 38 615-643. 
Ward, R. and Durrett, R. (2004). Subfunctionalization: How often does it occur? How 

long does it take? Theor. Popul. Biol. 66 93-100. 
Walsh, J. B. (1995). How often do duplicated genes evolve new functions? Genetics 139 

421-428. 

Watterson, G. a. (1983). On the time for gene silencing at duplicate loci. Genetics 105 
745-766. 

Wen, Z., Rupasinghe, S., Niu, G., Berenbaum, M. R. and Schuler, M. A. (2006). 

CYP6B1 and CYP6B3 of the Black Swallowtail [Papilio polyxenes): Adaptive evolution 

through subfunctionalization. Mol. Biol. Evol. 23 2434-2443. 
Wolfe, K. H. and Shields, D. C. (1997). Molecular evidence for an ancient duphcation 

of the entire yeast genome. Nature 387 708-713. 



Department of Mathematics 
Cornell University 
Ithaca, New York 14853 
USA 

E-MAIL: rtdl(amath.cornell.edu 



Department of Mathematics and Statistics 
Concordia University 
Montreal, Quebec 
Canada H3G IMS 

E-MAIL: lpopovic@mathstat. Concordia. ca 



